Three-dimensional adaptive evolution of gravitational waves in numerical relativity 
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Adaptive techniques are crucial for successful numerical modeling of gravitational waves from astro- 
physical sources such as coalescing compact binaries, since the radiation typically has wavelengths 
much larger than the scale of the sources. We have carried out an important step toward this goal, 
the evolution of weak gravitational waves using adaptive mesh refinement in the Einstein equations. 
The 2-level adaptive simulation is compared with unigrid runs at coarse and fine resolution, and is 
shown to track closely the features of the fine grid run. 
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A new era in general relativistic research is beginning 
with the inauguration of a worldwide network of grav- 
itational wave observatories including LIGO, VIRGO, 
GEO, and TAMA 0. The signals that these instruments 
will detect are expected to often arise in highly nonlinear, 
dynamical astrophysical events, such as the final merger 
of coalescing compact binaries. Fully general relativistic 
simulations of these events in three-dimensions will be 
an essential component of the successful detection and 
interpretation of gravitational wave signals. 

A key difficulty in calculating accurate waveforms from 
these simulations lies in the fact that the models en- 
compass both length and time scales that can vary by 
an order of magnitude or more. Adequate resolution 
is needed to model the strong-field dynamics near the 
sources. Also, the computational domain must be suf- 
ficiently large and well-resolved to model the resulting 
gravitational radiation that, for coalescing compact bi- 
naries, can have wavelengths ~ 10 — 100 times the size 
of the sources. Adaptive techniques that offer dynamic 
variable resolution are thus crucial for success. 

Adaptive mesh refinement (AMR) has undergone sig- 
nificant developments during the past two decades, and 
remains an active area of research and development. 
Within the context of general relativistic simulations, 
Choptuik first implemented the Berger-Oliger |^ AMR 
method to study critical phenomena in the collapse of 
massless scalar fields in spherical symmetry [^,|| ; see also 

. AMR techniques were applied to spherical black hole 
evolutions in Refs. 1^ and 0. In addition, the construc- 
tion of initial data for black hole collisions using AMR 
was studied in Ref. and fixed mesh refinement was 
employed for a short part of a binary black hole calcu- 
lation in Ref. ||]. Recently, Papadopoulos, Seidel and 
Wild |l^] carried out 3-D adaptive computation of grav- 
itational waves using a single model equation that de- 
scribes perturbations of a non-rotating black hole. 



We have carried out a 3-D simulation of pure gravita- 
tional waves in the Einstein equations using AMR. While 
our numerical relativity code does solve the full Einstein 
equations, we have chosen to evolve linearized Teukolsky 
waves to enable comparison with an analytic solution for 
this simulation. Overall, our 2-level adaptive simulation 
successfully reproduces the features of a uniform grid run 
at the finer resolution. 

Our numerical relativity code is based on the ADM 
code developed by the BBH Alliance and is writ- 
ten in Cartesian coordinates. We have implemented a 
conformal ADM formalism ||l^,^, and have incorpo- 
rated octant symmetry boundary conditions to minimize 
the number of gridpoints. Parallelism and AMR have 
been implemented using Paramesh , as discussed be- 
low. The time evolution is carried out using an iterative 
Crank-Nicholson method with two iterations [|l5| . Inter- 
polated Sommerfeld outgoing wave conditions are applied 
at the outer boundary of the grid |p^ , |l3| . 

Given the complexity of this system of equations and 
the computational techniques employed, analytic solu- 
tions provide an essential testbed for the resulting nu- 
merical code. To this end, we use the Teukolsky wave 
|l6| , which is a time-dependent solution to the linearized 
vacuum Einstein equations with geodesic slicing, a — 1 
and /3 = 0. It is based on a generating function F{t, r) — 
{A/uj'^){t ± r)exp(-(t ± r)'^/uj'^), where A « 1 IS an 
amplitude parameter and lu is the width of the gener- 
ating function, which is of the same order as the wave- 
length of the gravitational wave. To set initial data for 
our simulations, we use a time symmetric, even-parity 
L — 2, M — Teukolsky wave solution that contains a 
combination of ingoing and outgoing gravitational waves; 
see Refs. p6|Jl^]. We set = 1 and take A = 10"*^, 
so that the metric perturbation has an initial amplitude 
gij — 1 ~ ±10~^. Since t = is a moment of time sym- 
metry, we have Kij = initially. 
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This data is then evolved in our code with a = 1 and 
P = 0. The outgoing waves travel directly toward the 
edge of the grid. The initially ingoing waves first travel 
toward the origin, then reflect and move outward. As the 
overall signal travels away from the origin, it leaves flat 
space behind. Ref. gives explicit expressions for the 
metric components that we use as a first-order ana- 
lytic solution for comparison with our numerical results. 
All simulations presented in this paper were done with 
octant symmetry. 

We employ an AMR scheme that works on logically 
Cartesian, or structured, grids. The pioneers in this 
area are Berger and co-workers [^,^^, who covered 
the computational domain with structured grids and sub- 
grids that could overlap, take arbitrary shapes, be ro- 
tated with respect to the coordinate axes, and be merged 
with other sub-grids at the same level of refinement. Al- 
though this strategy is flexible and memory-efhcicnt, the 
resulting code can be very complex and difhcult to par- 
allelize, especially in 2D and 3D. Simplifications of this 
basic scheme were developed by Quirk ||2^ and later by 
DeZeeuw and Powell pl| ], who implemented refinement 
by bisecting the appropriate grid blocks in each coordi- 
nate direction and linking the hierarchy of sub-grids as 
the nodes of a data tree. Wild and Schutz have de- 

veloped a somewhat different approach using hierarchical 
linked-lists. 

Paramesh uses an AMR technique similar to that of 
DeZeeuw and Powell in which grid blocks are bi- 
sected in each coordinate direction when refinement is 
needed. The grid blocks all have the same logical struc- 
ture, with nxb zones in the x— direction, and similarly for 
nyb and nzb. Thus, in 3D refinement of a block yields 
8 child blocks, each having nxb x nyb x nzb zones but 
with zone sizes a factor of two smaller than in the par- 
ent block. Refinement can then continue on any or all of 
these child blocks, with the restriction that the grid spac- 
ing is not allowed to change by more than a factor of two, 
or one refinement level, at any location in the spatial do- 
main. Every grid block has a number of guard cell layers 
along each of its boundaries, to allow operations between 
neighboring blocks as well as between different refine- 
ment levels. These guard cells are filled with data from 
neighboring blocks or, if the block is along the edge of the 
computational domain, with appropriate outer boundary 
conditions. 

Paramesh handles the creation of grid blocks and 
builds and maintains the data structure (quad-tree in 
2D, oct-tree in 3D) that tracks the spatial relationships 
between blocks. It also handles all inter-block communi- 
cations. In addition, it keeps track of physical boundaries 
on which particular conditions are to be set, guarantee- 
ing that the child blocks inherit this information from 
the parent blocks. When executing on a parallel machine, 
Paramesh distributes the blocks among the available pro- 
cessors to achieve load balance, maximize block locality. 



and minimize inter-processor communication. 

We have implemented our numerical relativity code in 
the Paramesh framework. For the runs presented in this 
paper, we use blocks with nxb = nyb = nzb = nb and 
grid sizes Aa; = Ay — Az = h. We also use the same 
timestep, chosen to give stability on the finest grid, over 
the whole computational domain. The basic variables are 
located at the centers of the grid cells. We use second- 
order spatial finite differences and a single layer of guard 
cells. For the AMR runs, we use the restriction (transfer 
of data from fine to coarse grids) and prolongation (coarse 
to fine) operators built into Paramesh, which employ lin- 
ear interpolation. All the calculations presented in this 
paper were done using a T3E. 

A necessary preliminary step in using AMR for numer- 
ical simulations is to verify that the code works correctly 
when run using a single grid with constant h that covers 
the entire computational domain, i.e. for unigrid sim- 
ulations. Of particular importance is the convergence 
behavior of the code. Since we are using finite differ- 
ences that are second-order accurate in both space and 
time, we expect that the errors will decrease 0(h'^) as 

Using the Teukolsky wave initial data, we ran a set 
of three runs with increasingly finer grids: h = 0.25 (16'^ 
grid), h = 0.125 (32^ grid), and h = 0.0625 (64^ grid). In 
these runs, the outer boundaries of the grid were placed 
at Xmax = 2/max = -^max = 4. Calculation of various error 
norms demonstrated that, even for these relatively mod- 
est resolutions, the code has nearly second-order conver- 
gence for the propagation of the wave across the grid. 
However, as the waves move off the outer edge of the 
grid, lower amplitude reflected waves travel inward from 
the outer boundary. When the outer boundary is moved 
farther from the origin, and thus closer to the asymp- 
totic region in which the Sommerfeld outgoing wave con- 
ditions are applicable, the errors due to boundary effects 
decrease. 

To evolve gravitational waves using AMR, we have cho- 
sen to use a refinement criterion based on the first deriva- 
tive of the diagonal components of the conformal metric 
9ij — ^~^'^9iji where e'*'^ = dei{gijY/^ . Thus, in each 
zone we calculate the quantity t — ^(^§^ + §^ + §7)1 
where / = cjij — 1 for i = j. The computational do- 
main is covered with blocks having nb^ physical zones. 
A block is marked for refinement if r > for at least 
one zone in that block, and for derefinement if r < Tdo for 
all zones in that block, where Tre > T^e. Note that both 
refinement and derefinement in Paramesh are restricted 
by the requirement that the grid can change by only one 
refinement level at block boundaries. 

We concentrate on 2-level AMR, with level 1 indicating 
the coarse grid and level 2 the fine grid. For the simu- 
lation presented here, we use grid blocks with nb — A 
and set the outer boundary of the computational domain 
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at Xinax = Umax = ^^max = 7. After some experimenta- 
tion, we chose to trigger refinement by r^e = 1.75, and 
derefinement by r^e = 0.875. 

To set up the initial data, we first covered the entire 
computational domain with level 1 blocks, and used the 
Teukolsky solution to give gij and (f> . The quantity r 
was then computed and the refinement criterion applied, 
yielding blocks at level 2. The Teukolsky solution was 
again applied to these refined blocks to set up gij and 4>. 
However, when the first derivatives of cjij were taken nu- 
merically to compute the conformal connection functions 
r* discontinuities in F* appeared at the boundaries 
between level 1 and level 2 blocks. When this data was 
then evolved forward in time, these initial discontinuities 
produced errors in other variables. To circumvent these 
problems, we calculated on a level 2 grid covering the 
whole domain, and then used the restriction operator to 
give F' on the level 1 blocks. This produced a data set 
without initial discontinuities. 

Figure |l| shows the metric component g^z hi the equa- 
torial plane for an AMR run with resolution h = 0.25 
on level 1 and h = 0.125 on level 2 at time t = 5.06. 
Initially, only the region near the origin is covered by 
the fine grid. This refined region moves outward as the 
wave propagates toward the edge of the grid. As seen 
m Fig. the region near the origin is derefined as the 
solution approaches flat space in that region. 

It is interesting to compare this AMR solution with 
unigrid runs using the same computational domain. Fig- 
ure |2| shows the quantity gzz — 1 along the x— axis at 
t = 5.06, with the data from the AMR run shown as 
filled squares. In addition, a short dashed line shows 
the analytic solution, a solid line the unigrid result with 
h = 0.125 (56^ grid, same resolution as the level 2 AMR 
grid) and a long dashed line the unigrid result with 
h = 0.25 (28'^ grid, same resolution as the level 1 AMR 
grid). Notice that the AMR solution closely tracks the 
higher resolution unigrid result. When we take the Li, 
L2, and Loo error norms of gzz over the entire computa- 
tional volume, we find that the AMR results closely track 
those of the higher resolution unigrid run. 

Initially, less than 10% of the computational volume 
was covered by the fine grid. This fine grid volume frac- 
tion increases as the wave propagates outward, as shown 
in Fig. ^. Although derefinement in the equatorial plane 
first occurs around the origin at t ~ 5, the total number 
of fine grid cells continues to increase as the expanding 
wavefronts propagate outward. The main peak of the 
wave (located at 2; ~ 4.5 in Fig.|^) encounters the outer 
boundary at i 7.5. Around this time, Fig. |^ shows 
that the fine grid volume fraction begins decreasing as 
the derefined region behind the wave grows. By i ~ 9, 
most of the Teukolsky wave has left the grid. However, 
refiected waves traveling inward from the outer boundary 
cause continued refinement and derefinement, preventing 
the fine grid volume fraction from dropping to zero. 



In this calculation, the savings in computational cost 
of the AMR run over a unigrid run are modest, even if 
we neglect the problems with the boundary conditions, 
since the gravitational wavelength A is comparable to the 
dimension L of the computational domain. In a more re- 
alistic simulation, we would expect A <C i, resulting in 
significant computational savings. For a full simulation 
of binary coalescence, one will need to handle multiple 
scales to evolve both the dynamics of the sources and 
the propagation of the waves that develop. While the 
sources are expected to require high resolution, particu- 
larly near regions of shocks or shear layers, the waves can 
likely be handled with lower resolution. Thus, although 
the resolutions used here are modest compared to those in 
previous numerical studies of Teukolsky waves, e.g. 
they are typical of the resolutions likely to used for grav- 
itational waves in full simulations of binary coalescence 
in the next few years. 

While this AMR calculation consitutes an important 
step toward the ultimate goal of realistic models of gravi- 
tational wave sources, more work remains to be done. For 
example, even though AMR can be used to alleviate the 
problem of reflected waves by moving the outer bound- 
ary even farther from the origin, better outer boundary 
conditions are still needed |^^. In addition, it will be 
interesting to test different refinement criteria, includ- 
ing those based on truncation error estimation [p|,p^. 
Finally, we observed that derefinement introduces some 
high frequency noise, principally in the extrinsic curva- 
ture variables. While this did not pose a significant prob- 
lem for this 2-level run, it can become important when 
more refinement levels are used. We are experimenting 
with various ways of removing this noise by using filter- 
ing or adding dissipation [ pSp^ , and will report on 
this work elsewhere. 
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FIG. 2. The quantity Qzz — 1 is shown along the a;— axis at 
time t = 5.06. The filled squares show the data for the 2-level 
AMR run, the short dashed line the analytic result, the solid 
line the unigrid result with h = 0.0625, and the long dashed 
line the unigrid result with h = 0.125. 
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